####################################################################################
###
###   Immigration, Public Housing and Support for the French National Front
###   Gloria Gennaro
###
###   Paper Table A3, A4
###
####################################################################################


library(multiwayvcov)

################################################################################
# Set Up
################################################################################

# Load Working Directories
source(paste0(wd_main, '/2_code/00_working_directories.R'))

# Load data
df = read.csv(paste0(wd_data, '/dataset_newspaper_analysis.csv'))

################################################################################
# Auxiliary Functions
################################################################################

LinearCombTest = function (lmObject, vars, .vcov = NULL) {
  if (is.null(.vcov)) .vcov = vcov(lmObject)
  beta = coef(lmObject)
  sumvars = sum(beta[vars])
  se = sum(.vcov[vars, vars]) ^ 0.5
  tscore = sumvars / se
  pvalue = 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
  matrix(c(sumvars, se, tscore, pvalue), nrow = 1L,
         dimnames = list(paste0(paste0(vars, collapse = " + "), " = 0"),
                         c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
}

lincomCoeff = function(model){
  out = data.frame()
  vcv = cluster.vcov(model, cluster = df$dep)
  out = rbind(out, c(LinearCombTest(model, c("x"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x", "x:factor(imm_quant_99)2"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x", "x:factor(imm_quant_99)3"), vcv), nobs(model)))
  colnames(out) =  c('coef', 'se', 'tval', 'pval', 'nobs')
  return(out)
}

marignaleffects = function(model){
  out = data.frame()
  vcv = cluster.vcov(model, cluster = df$dep)
  out = rbind(out, c(LinearCombTest(model, c("x"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x:factor(imm_quant_99)2"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x:factor(imm_quant_99)3"), vcv), nobs(model)))
  colnames(out) =  c('coef', 'se', 'tval', 'pval', 'nobs')
  return(out)
  
}



################################################################################
# Estimation and plot functions
################################################################################

risultati = function(name){
  
  pars = as.list(match.call()[-1])
  
  # Select variables
  df$y = scale(df[,as.character(pars$name)]/df$counts_newspapers)
  df$x =  scale(df$sruT)
  
  # Main analysis
  container_lincom = data.frame()
  container_marg = data.frame()
  
  formula = 'y ~ x*factor(imm_quant_99) + factor(dep) + factor(year)'
  controls = c('', '+ factor(reg)*year', '+ factor(reg)*year + log(pop_dep)',
               '+ factor(reg)*year + log(pop_dep) + I(num_res_own/num_log) + I(num_empty_log/num_log)')

  for (ctr in controls){
    model = lm(as.formula(paste0(formula, ctr)), data=df)
    container_lincom = rbind(container_lincom, lincomCoeff(model))
    container_marg = rbind(container_marg, marignaleffects(model))
  }
  
  # Linear Combinations 
  x = scale(df$sru_weight)
  model = lm(y ~ x*factor(imm_quant_99) + factor(dep) + factor(year) + factor(reg)*year, data=df) # %>% summary()
  container_lincom = rbind(container_lincom, lincomCoeff(model))
  container_marg = rbind(container_marg, marignaleffects(model))
  
  # Clean container
  clean_container = function(container){
  colnames(container) = c('coef', 'se', 'tval', 'pval', 'nobs')
  container$group = rep(c('Low\nImmigration', 'Medium\nImmigration', 'High\nImmigration'), 5)
  container$estimation = rep(c('No Controls', '+ Regional Trends', '+ Population', '+ Housing market', '+ Weighted'),each=3)
  container$group.f = factor(container$group, levels = rev(unique(container$group)))
  container$estimation.f = factor(container$estimation, levels = unique(container$estimation))
  container_stored = container
  }
  
  
  container_lincom = clean_container(container_lincom)
  container_marg = clean_container(container_marg)
  output = list('lincom' = container_lincom, 'margins' = container_marg)
  return(output)

}


cleancontainer = function(cont){
  container = cont
  nobs = container$nobs
  container$se = paste0('(', as.character(round(container$se, 3)), ')')
  container$coef = as.character(round(container$coef, 3))
  container$coef = ifelse(container$pval>0.1, container$coef,
                          ifelse((container$pval<0.1 & container$pval>0.05), paste0(container$coef, '*'),
                                 ifelse((container$pval<0.5 & container$pval>0.01), paste0(container$coef, '**'),
                                        paste0(container$coef, '***'))))
  
  # Build latext table
  container = container[c("coef", "se", 'group', 'estimation')]
  
  coeffs = paste0(container[which(container$estimation == "No Controls" ), c("coef")], ' & ',
                  container[which(container$estimation =="+ Regional Trends" ), c("coef")], ' & ',
                  container[which(container$estimation =="+ Population" ), c("coef")], ' & ',
                  container[which(container$estimation =="+ Housing market" ), c("coef")], ' & ',
                  container[which(container$estimation =="+ Weighted" ), c("coef")])
  
  ses = paste0(container[which(container$estimation == "No Controls" ), c("se")], ' & ',
               container[which(container$estimation =="+ Regional Trends" ), c("se")], ' & ',
               container[which(container$estimation =="+ Population" ), c("se")], ' & ',
               container[which(container$estimation =="+ Housing market" ), c("se")], ' & ',
               container[which(container$estimation =="+ Weighted" ), c("se")])
  
  nobs = paste(nobs, collapse = ' & ')
  
  output = paste0(paste(paste0(c('Low Immigration', 'Medium Immigration', 'High Immigration'), '&', coeffs, '\\', ses), collapse='\\'),
                  '\\', nobs)
  return(output)
}


################################################################################
# Estimates
################################################################################

out1 = risultati(score)
out1_neg = risultati(score_neg)


################################################################################
# Table A3
################################################################################

print(cleancontainer(out1$margins)) # panel 1
print(cleancontainer(out1$lincom)) # panel 2

################################################################################
# Table A4
################################################################################

print(cleancontainer(out1_neg$margins)) # panel 1
print(cleancontainer(out1_neg$lincom)) # panel 2
